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Abstract 

We use single-cluster Monte Carlo simulations to study the role of 
topological defects in the three-dimensional classical Heisenberg model 
on simple cubic lattices of size up to 80'^. By applying reweighting 
techniques to time series generated in the vicinity of the approximate 
infinite volume transition point Kc, we obtain clear evidence that the 
temperature derivative of the average defect density d{n)/dT behaves 
qualitatively like the specific heat, i.e., both observables are finite in 
the infinite volume limit. This is in contrast to results by Lau and 
Dasgupta [Phys. Rev. B39 (1989) 7212] who extrapolated a divergent 
behavior of d{n) / dT at Kc from simulations on lattices of size up to 16^ . 
We obtain weak evidence that d{n)/dT scales with the same critical 
exponent as the specific heat. As a byproduct of our simulations, we 
obtain a very accurate estimate for the ratio a/u oi the specific-heat 
exponent with the correlation-length exponent from a finite-size scaling 
analysis of the energy. 
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1 Introduction 



It is well known that topological defects can play an important role in phase 

transitions.^'^ Extensively studied examples of systems with pointlike de- 
fects are the two-dimensional (2D) XY modeF and defect models for 2D 
melting.^'^ Recently Lau and Dasgupta (LD)^ have used Monte Carlo (MC) 
simulations to study the role of topological defects in the three-dimensional 
(3D) classical Heisenberg model, where the defects are also point-like objects 
with a binding energy that increases linearly with the separation.^ Moti- 
vated by the importance of vortex points in the 2D XY model, LD tried to 
set up a similar pictorial description of the phase transition in the 3D Heisen- 
berg model. Analyzing their simulations on simple cubic (sc) lattices of size 
V — with L = 8, 12 and 16, LD claimed that the temperature deriva- 
tive of the average defect density, (n), diverges at the critical temperature 
Tc like d{n)/dT ~ t~^, t = \T — Tc\/Tc, with an exponent ip ~ 0.65. They 
further speculated that ip = 1 — /3, where /9 ~ 0.36 is the critical exponent of 
the magnetization, and then argued that (n) should behave like a "disorder" 
parameter. 

At first sight, the existence of such a strong divergence of d{n)/dT seems 
unlikely, because the definition of defects is quasi-local. It is therefore more 
likely^ that (n) should qualitatively behave like the energy and d{n)/dT like 
the specific heat, which is a finite quantity for the 3D Heisenberg model. 

Using standard finite-size scaling (FSS) arguments we hence expect to see 
on finite lattices either 

d{n)/dT ^ L'^/^'fix) (1) 

or, if the second argument holds true, 

d{n) /dT = const + L°/'^^(x) , (2) 

where u ~ 0.7 and a ~ —0.1 are the correlation length and specific heat expo- 
nents, respectively, x = tL^/^ is the finite-size scaling variable, and f{x),g{x) 
are scaling functions. At fixed x, Ansatz (1) predicts an approximate linear 
divergence in L, d{n)/dT ~ L'^^'^ L. On the other hand, because a is nega- 
tive for the 3D Heisenberg model, Ansatz (2) predicts a constant asymptotic 
value, d{n)/dT ~ const. For sufficiently large L a clear distinction between 
(1) and (2) should hence be observable. 

To decide between the two alternatives we have performed MC simula- 
tions on large lattices of size up to L = 80, employing the single-cluster 
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update algorithm^ and reweighting techniques.^ As a result we find unam- 
biguous support for tlie second alternative, d{n)/dT ~ const, as L — > oo. 

As a byproduct of our simulations we numerically extrapolate a/v from 
a FSS of the energy in close vicinity of the critical coupling Kc with a much 
higher accuracy than obtained in recent high-precision MC studies.^°'^^ 



2 The simulation 

The partition function of the Heisenberg model is given by 



(Elj 



47r 



eM-KE)- E = Y,{l-s,-s,). (3) 



where K = J/ ksT is the (reduced) inverse temperature, Si are three-dimensional 
unit vectors at the sites i of a sc lattice, and denotes nearest-neighbor 
pairs. Using the single-cluster update algorithm^ we have simulated the par- 
tition function (3) for lattices of size V = with L=8,12, 16,20, 24, 32, 40, 
48, 56, 64, 72, 80 and periodic boundary conditions. Our main emphasis was 
on the defect density n = Z]q'^n|q|, where ni,n2, . . . are defect densities of 

charge q — ±1, ±2, To locate these charges we followed the definition of 

Berg and Liischer^^ according to which the charge qi* at the dual lattice site 
i* is given by 

?^* = ^E^- (4) 



cos(iA) = l + s^.s, + s,.ss + ss-s, 

2 J2(l + Si • S2){1 + S2 ■ 53)(1 + 53 • Si) 



with 



and the sign of Ai is determined by signAj = sign (si • (s2 x ss)). For a sc 
lattice the distinction between i and i* is inessential, since the difference 
is only a uniform translation along the space diagonal. The sum J2i=i in 
(4) refers to the 12 triangles that can be formed on the faces of the cube 
enclosing qi*. In (5) the spins at the corners of these triangles are numbered 
in a counter-clockwise sense relative to the outward pointing normal. For the 
orientation of the diagonals on each face of the cube we used the convention 
that they run from the lattice point i to i + ei + ej, and from i + ei to 
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i + Si + ej + ek, where e.;, i = 1, 2, 3 denote unit-vectors in the three coordinate 
directions. This choice is obviously not unique, but we have checked that 
other conventions give on the average the same charges within the error 
bars. From the definition of qi* it is clear, that a trivial upper bound on the 
magnitude of the lattice topological charge is g^. < 5. In our runs the highest 
topological charge observed was three, which occured on the order of 10"''' 
per site and measurement, see Table 1. The likelihood of the appearence of 
the higher charges was probably too small for them to occur during our run 
times. 

All runs were performed close to the approximate infinite volume tran- 
sition point Kc = 0.6930, as determined in recent MC studies^°'^^'^^ of this 
model. Since we wanted to have reference data we performed our simulation 
at the same coupling Kq — 0.6929 as in our study of ref.,^^ which is close 
enough to Kg. to allow safe reweighting of our data to this value of Kc- Be- 
cause the computation of qi* is quite complex and thus time consuming, we 
have performed many cluster update steps between measurements, adjusted 
in such a way that the (integrated) autocorrelation time of the charge density 
measurements is around t„ ~ 1 — 2. Since it turned out that the (integrated) 
autocorrelation times t„ and of the charge density and the magnetic sus- 
ceptibility are roughly equal, we were able to guess the required measurement 
interval by extrapolating our previous results for t^^^ to larger lattice sizes 
L. The measurement statistics are given in Table 1. While the statistics is 
comparable to that of our previous studies, and much better than that of 
LD, we note, that our investigated lattices have much larger linear size up 
to L=80 as compared to L=48 in our previous work, and compared to the 
largest size L = 16 of LD. For each run we recorded the time series of the 
energy density e = E/V, the magnetization density m = \J2iSi\/V , and the 
charge densities n\q\. The resulting averages (e), (n), and {n\g\) can be found 
in Table 1. 

To compute the specific heat C = d{e)/dT, the thermal expansion co- 
efficient Cq = Td{n) /dT, and the topological susceptibility Xq = d{n)/diJ,, 
where //. is the "field" in a fugacity term //, J2i* Qi* which one can imagine 
adding to the energy in (3), and n is defined as n = rii -|- 4^2 + 9^3 -|- . . ., we 
used the relations 

C = VK\{e')-{e)') = VK'{e;e), (6) 
Cq = VK{{en)-{e){n)) = VK{e;n), (7) 
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X, = V{{n')-{nf) = V{n;n). 



(8) 



To obtain results for the various observables O dX K values in an inter- 
val around the simulation point Kq = 0.6929, we applied the reweighting 
method.^ Since we recorded the time series this amounts to computing 



(Ce 



-AKE\ 



(^)u= 7:-A..., (9) 



-AKE 



\Ko 



with AK = K — Kq. To obtain errors we devided each run into 20 blocks 
and used standard Jackknife errors. 

The results for the quantities in (6) - (8) at Kc = 0.6930 arc collected in 
Table 2. Also given are the eigenvalues Ai, A2 of the 2x2 covariance matrix 
M of e and n with elements Mn = VK'^{e;e), M12 = M21 = VK{e;n), 
M22 = V{n;n). 



3 Results 

Applying (9) we have determined the temperature dependence of the quan- 
tities in (6)- (8). For small lattices, Cg has its peak location at temperatures 
larger than Tc, in contrast to C, which peaks at temperature values smaller 
than Tc. With increasing lattice size, however, we observe a strong correla- 
tion between Cg and C, that is, both quantities develop a smooth peak at 
roughly the same temperature (T < Tc), see Fig. 1. In contrast to C the 
peak locations of Cg scale non-monotonically, with a crossover at L a; 20. 

We focussed first on the scaling behavior of Cg at our previous estimate of 
the critical coupling Kc — 0.6930, obtained from the crossings of the Binder 
parameter U — 1 — {m'^) /3{m,^)^}^ Our new data for U on the large lattices 
confirmed the constancy of our previous result of U* = 0.6217(8) and hence 
our estimate for Kc- We checked first a scaling Ansatz for Cg of the form 

C, = - aoL">, (10) 

where Cg"^^ is a regular background term which is assumed to be independent 
of lattice size.-*^^ Note that this Ansatz covers both scahng hypotheses (1) 
and (2). The resulting fit shown in Fig. 2(a) yields a'/u = —0.401(61), 
C'^^s = 1.50(8), and ao = 1.82(6), with a quality factor^^ Q = 0.30. The good 
quality of the fit basically rules out the divergence predicted by the Ansatz 
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(1) of LD, and strongly favours (2), which predicts a finite asymptotic value 
for Cq. Only if one assumes that the FSS behavior sets in at extremely large 
lattice sizes, one could still attain an assertion of the form (1), but with the 
consequence of an extremely small exponent 'ip. We also tried to reproduce 
the exponent ~ 0.65 of LD, by selecting only their lattices sizes, and 
fitting a straight line to our first 3 data points. But even then we obtain a 
much smaller value of ■(/'/z/ ~ 0.36(3), leading io ip k, 0.25(3). We think, this 
discrepancy to the result of LD is partly due to our higher statistics and partly 
due to the fact that we obtained Cq through a thermodynamic derivation, 
which normally gives better results than the numerical differentiation used 
by LD. 

Because the Ansatz (2), which was based on the assumption that (n) 
should behave hke the energy, fits so well, one can ask, if a' is equal to the 
specific-heat exponent a. Using our earlier MC result^^ of u — 0.704(6), we 
get a value of a' = —0.282(46), which does, on the first glance, not strongly 
support this conjecture. The best field theoretical estimates arc u = 0.705(3), 
a = —0.115(9), and ajv = —0.163(12) (resummed perturbation series^^), 
while our earlier MC study^^ yielded v = 0.704(6), a = -0.112(18), and 
ajv — —0.159(24). However, the accuracy of the values of a is somewhat 
misleading, because they were obtained from hyperscaling, a = 2 — 3z/. The 
directly measured values have much larger error bars, for example ajv — 
-0.30(6)^° and ajv = -0.33(22). 

To compare a' directly with the measured specific-heat exponent of the 
present MC simulation, we fitted C to 

C = C'^s - 6oL"/^ (11) 

with a constant background term C^*^^.^^ The resulting fit in Fig. 2(b) yields 
a/u = -0.225(80), = 4.8(7), and bo = 4.1(5) with Q = 0.55, leading 
to a = —0.158(59). These values are in very good agreement with the hy- 
perscaling prediction, but noteworthy is also the tendency for the values to 
come out too large. 

Of course, a fit of a divergent quantity, like the first derivative dCg/dT, 
for example, is in principle numerically much easier to handle. We tried to do 
this for dCq/dT and dC/dT at K^, and observed the expected divergent scal- 
ing behavior, but unfortunately the statistical errors of the third cumulants 
involved turned out to be much too large to allow for meaningful fits. 
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Other estimates for a and a' can be obtained by means of fits of (e) and 
(n), which again look quahtatively very ahke. According to (11), we should 
have on periodic lattices a scaling behavior of the energy density (e) of the 
form^^'^^ 

(e) = {ef^^ - d^L^^'-^^l^ (12) 

and because of (10) the topological charge density (n) should then accord- 
ingly scale like 

(n) = {nY^^ - coL^"'-^)/". (13) 

Fits of these quantities at = 0.6930, shown in Fig. 3, yield (a' — l)/z^ = 
-1.547(15), (n)^*^g = 0.1074(1), and Cq = 0.42(2), with Q = 0.30, and {a - 
= -1.586(19), (e)''^s = 2.0106(1), and do = 1.68(8), with Q = 0.25. 
This results in a'/iy = -0.127(27), a' = -0.089(20), a/u = -0.166(31), and 
a = —0.117(23). The results for a and a/v are in excellent agreement with 
the hyperscaling prediction, and have not been directly measured before with 
such a high precision. We attribute this to our large lattice sizes used, but 
also to the fact that we used fits of (e) instead of C. The results for a' and 
a'/u are now lower than those obtained in (10), but now they are almost 
consistent with the values for a and a/u. Still it is a little bit puzzling that 
both estimates for the exponent a'/u obtained from the fits (10) and (13) 
do not agree in their respective error range. We attribute this partly to the 
unknown FSS behavior of the regular background term (n)""*^^, and partly to 
the fact that the statistical errors of the three parameter fits should be taken 
with great care. 

We further looked at the scaling behavior of Xg, defined in eq.(8). A first 
look at the plots suggests to try again a scaling Ansatz of the form 

X, = - eoi^""/^ (14) 

From a three-parameter fit we obtain a"/]^ = —0.554(57), xl^^ = 0.67(2), 
and Co = 0.95(6) with Q = 0.41, leading to a" = -0.390(44). This time it 
seems already very unhkely, that a" is equal to the specific-heat exponent. 
However, if one discards the two lowest L values from the fit, one observes 
a clear trend towards a lower Q;"-value, but with the drawback of increased 
error bars and no improvement in x^/dof (per degree of freedom). 

We also checked in all other fits if there were corrections to FSS, by 
discarding successively the data points for L = 8 and L = 12. We observed 
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in all quantities a trend to the value ol a/v predicted by hyperscaling, but 
at the price of much larger error bars. Also the x^/dof did not improve. We 
further checked for confluent corrections/^ by including a term of the form 
OiL"'^, with to — ^/v fixed at the literature value ^.1%}"^ But again the fits 
were too unstable to give conclusive results. 

We also tested if our results depended strongly on the choice of K^., by 
repeating the fits of all quantities at K(. ± 0.0002. The resulting parameters 
were always consistent with the values at Kc in the one-o" range. 

To get a clearer picture we further looked at the scaling behavior of the 
eigenvalues of the covariance matrix M of e and n, defined by equations (6) 
- (8), which give two uncorrelated observables Ai and A2. Again we used a 
scaling Ansatz of the form 

A, = Ar^-a,W^ (15) 

As results we obtain a^/v = -0.273(73), \'^^ = 5.1(5), and Oi = 4.7(2), 
with Q = 0.49 and a2/v = -1.45(42), X'^^ = 0.1307(8), and 02 = 0.2(2), 
with Q = 0.60, leading to ai = -0.192(54) and ^2 = -1.02(31). This 
suggests ai ^ a and 02 ~ a — 1. The quality of the fits can be inspected 
in Fig. 4. The existence of an uncorrelated observable which scales with an 
exponent different from a suggests that there is possibly a new scaling field, 
and that Cq and Xg see remnants of this scaling field in their FSS behavior. 
This becomes particularly clear if one remembers that Ai + A2 = C + Xg- 
Therefore at least Xg should see something of the exponent 0:2 ~ a — 1. 
Another alternative would be, that Cg and Xg do not scale directly with a/u 
but with some rational multiple a/v. As long as there is no satisfactory 
theory of the scaling of topological quantities, however, one cannot decide 
between these alternatives. 

4 Concluding remarks 

We have shown that in the three-dimensional classical Heisenberg model 
the topological defect density (n) and its temperature derivative Cq behave 
qualitatively like the energy (e) and its temperature derivative C. Especially, 
we can reject the conjecture of LD that Cq diverges with a new critical 
exponent ■0. 
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Rather, our simulations indicate that Cq behaves also quantitatively like 
the specific heat, i.e. scales like (2). We obtain weak evidence that asymp- 
totically for large L the scaling of Cq is governed by the specific-heat critical 
exponent a. Still, it cannot be ruled out that the scaling of Cq involves also 
a new exponent belonging to a new scaling field. For the topological suscep- 
tibility Xq wc find that it also remains finite, and that it can be fitted with an 
Ansatz of the form (2) as well, but that its scaling exponent must be some 
multiple of a, be the representative of a new scaling field or be a mixture of 
both. Our fits of the eigenvalues Aj of the covariance matrix seem to indicate 
that Cg and mixture of a part which scales with a, and a part which 

scales according to a — 1. 

Finally, the present fits of the specific heat at Kc yielded a value of a of 
better accuracy and in better agreement with the hyperscaling value than 
fits of the specific-heat maxima as used in previous works, which we 
attribute to our large lattice sizes, the larger number of available data points, 
and to the fact, that our data and fit was done cxtrcmly close to the critical 
temperature. Moreover, by fitting the energy at Kc to eq. (12), we obtained 
an estimate for a/u with a precision so far unpreceeded by direct numerical 
MC simulations and in accuracy comparable to hyperscaling predictions. 
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Tables 



L 


iVo 


TV 

meas 




(e) 


(n) X 10 


(ni) X 10 


(n2) X 10^ 


(na) X 10^ 


8 


17 


50178 


1.2 


1.9487(9) 


0.9054(18) 


.88090 


6.13 


7.78 


12 


20 


159575 


1.6 


1.9786(3) 


0.9845(7) 


.95689 


6.90 


5.80 


16 


40 


64368 


1.2 


1.9905(3) 


1.0170(9) 


.98838 


7.16 


7.21 


20 


50 


27670 


1.3 


1.9968(3) 


1.0347(6) 


1.0053 


7.35 


5.42 


24 


50 


20000 


1.5 


1.9998(2) 


1.0431(6) 


1.0136 


7.38 


9.04 


32 


68 


25403 


1.5 


2.0045(2) 


1.0561(4) 


1.0260 


7.52 


6.61 


40 


74 


21765 


1.9 


2.0063(1) 


1.0617(3) 


1.0314 


7.59 


8.04 


48 


93 


21005 


1.9 


2.0074(1) 


1.0646(3) 


1.0342 


7.61 


7.83 


56 


136 


23795 


1.6 


2.0084(1) 


1.0674(2) 


1.0369 


7.62 


6.56 


64 


200 


26439 


1.4 


2.0090(1) 


1.0691(1) 


1.0385 


7.64 


7.18 


72 


150 


20000 


1.8 


2.0093(1) 


1.0701(2) 


1.0395 


7.65 


7.10 


80 


200 


25431 


1.7 


2.00962(4) 


1.0709(1) 


1.0403 


7.66 


6.77 



Table 1: Measurement statistics at the simulation point — 0.6929: L is 
the linear lattice size, Nq is the number of cluster steps between measure- 
ments, A'^meas stauds for the number of measurements, t„ is the integrated 
autocorrelation time of the charge density, (e) is the energy density, {n\q\) 
are the observed densities of dual cells with charge \q\, and the total defect 
density {n) is defined as {n) — {ni) + 4(^2) + 9(^3) + 
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L 


C 


c. 


Xq 


Ai 


A2 


8 


2.177(18) 


0.712( 6) 


0.369(3) 


2.423(19) 


0.1225( 8) 


12 


2.407(12) 


0.833( 5) 


0.430(2) 


2.721(14) 


0.1261( 6) 


16 


2.562(27) 


0.908(10) 


0.467(4) 


2.900(30) 


0.1283( 9) 


20 


2.651(36) 


0.951(14) 


0.486(6) 


3.009(40) 


0.1277( 7) 


24 


2.752(40) 


0.993(17) 


0.506(8) 


3.128(47) 


0.1297(12) 


32 


2.832(38) 


1.027(14) 


0.521(6) 


3.223(43) 


0.1308(12) 


40 


2.970(45) 


1.081(16) 


0.542(7) 


3.382(51) 


0.1300( 9) 


48 


2.977(44) 


1.097(18) 


0.553(8) 


3.400(51) 


0.1304(13) 


56 


3.142(35) 


1.157(12) 


0.576(5) 


3.587(38) 


0.1309(13) 


64 


3.181(30) 


1.173(10) 


0.579(5) 


3.631(33) 


0.1286(10) 


72 


3.142(55) 


1.155(19) 


0.574(7) 


3.585(62) 


0.1315(14) 


80 


3.182(49) 


1.169(16) 


0.579(6) 


3.630(53) 


0.1307(12) 



Table 2: Results for the matrix elements of the covariance matrix V{Oi] Oj), 
with Oi = Xe, O2 = n at X = 0.6930(Ri Kc). Also included are the 
eigenvalues Ai, A2 of the covariance matrix. 



12 



Figure Headings 



Fig. 1: (a) Cq = Td{n)/dT, and (b) the specific heat C versus K for lattices 
of size L = 12, 40, 72. The values were obtained by reweighting of the 
runs at Kg = 0.6929. 

Fig. 2: (a) Cg, (b) C, and (c) at = 0.6930(ai K^) as function of the 
lattice size L. The solid lines show the best non-linear three-parameter 
fits to the data. 

Fig. 3: (a) (n), and (b) (e) at K = 0.6930(Ri Kc) as function of L, together 
with the best non-linear three-parameter fits to the data. 

Fig. 4: (a) Ai and (b) X2 Sit K — 0.6930(Ri Kc) as function of L. Included 
are also the best non-linear three-parameter fits to the data. 
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